%% Function to compute the equilibrium in rural areas
function f = RuralEquilibriumFunction(pp,pop_R, XS,XA,alphaA,alphaM,alphaS,cAbar,cSbar, Xn,Xr,XrXr,XnXn,f_XrXn_R);

   
%% Unpackaging prices   
    p_SR    =  pp;
    wu=XA;

%% Constructing measures of skills distributions in Urban Regions (U= C + F)    
   
    Xn_f_XrXn_R   =  XnXn.*f_XrXn_R;
    Xr_f_XrXn_R   =  XrXr.*f_XrXn_R;

% Implied wages and labor market earnings
    wrR                     =    p_SR*XS;
    ybar_R                  =    p_SR*cSbar*(1-alphaS)/alphaS+cAbar;       % income threshold for the consumption of Services (presumes zero housing costs in rural areas)
% Labor Markets
    yR                      =    max(wu,wrR*XrXr);
    IOCC_r_R                =    zeros(size(XrXr));  % indicator, occupation choices
    IOCC_r_R(wu<wrR*XrXr)   =    1;
    ICS_R                   =    zeros(size(XrXr));  % Indicator, positive 
    ICS_R(yR>=ybar_R)       =    1;
    LrR                     =    trapz( Xr,trapz(Xn, Xr_f_XrXn_R .*IOCC_r_R));  % effective labor in routine tasks, rural areas
    ErR                     =    trapz( Xr,trapz(Xn,    f_XrXn_R .*IOCC_r_R));  % employment in routine tasks, rural areas
    QA                      =    XA*(pop_R-ErR);
    QS_R                    =    XS*LrR;
% Consumption of Services
    yTop_R                  =    trapz( Xr,trapz(Xn, yR.*f_XrXn_R.*ICS_R));
    fTop_R                  =    trapz( Xr,trapz(Xn,     f_XrXn_R.*ICS_R));
% Equilibrium Pricing Condition
    PP_SR1                  =    alphaS*(yTop_R-cAbar*fTop_R)/(QS_R +(1-alphaS)*cSbar*fTop_R);

%% Distance (Euclidean)

    f = (PP_SR1-p_SR)^4;  
 
 
